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We  use  the  recently  developed  reactive  force  held  ReaxFF  with  molecular  dynamics  to  study 
thermal  induced  chemistry  in  RDX  [cychc-[CH2N(N02)]3]  at  various  temperatures  and  densities. 
We  hnd  that  the  time  evolution  of  the  potential  energy  can  be  described  reasonably  well  with  a 
single  exponential  function  from  which  we  obtain  an  overall  characteristic  time  of  decomposition 
that  increases  with  decreasing  density  and  shows  an  Arrhenius  temperature  dependence.  These 
characteristic  timescales  are  in  reasonable  quantitative  agreement  with  experimental  measurements 
in  a  similar  energetic  material,  HMX  [cychc-[CH2N(N02)]4] .  Our  simulations  show  that  the 
equilibrium  population  of  CO  and  CO2  (as  well  as  their  time  evolution)  depend  strongly  of  density: 
at  low  density  almost  all  carbon  atoms  form  CO  molecules;  as  the  density  increases  larger 
aggregates  of  carbon  appear  leading  to  a  C  dehcient  gas  phase  and  the  appearance  of  CO2 
molecules.  The  equilibrium  populations  of  N2  and  H2O  are  more  insensitive  with  respect  to  density 
and  form  in  the  early  stages  of  the  decomposition  process  with  similar  timescales.  ©  2005 
American  Institute  of  Physics.  [DOI:  10.1063/1.1831277] 


I.  INTRODUCTION 

A  molecular  level  understanding  of  condensed-matter 
chemistry  is  a  central  problem  in  many  areas  of  chemistry 
and  materials  science.  Despite  the  enormous  experimental 
and  theoretical  progress  in  recent  years  a  wide  variety  of 
phenomena  remains  uncharacterized.  Particularly  challeng¬ 
ing  is  understanding  the  decomposition  and  subsequent  reac¬ 
tions  of  high  energy  (HE)  materials  due  to  the  extreme  con¬ 
ditions  of  temperature  and  pressure  involved  as  well  as  their 
complex  chemistry.  Recent  advances  in  experimental  tech¬ 
niques  such  as  ultra-fast  spectroscopy^  and  interferometry^ 
are  beginning  to  provide  molecular-level  information  about 
the  decomposition  of  HE  materials.  Nevertheless,  molecular 
dynamics  (MD)  remains  the  only  technique  capable  of  pro¬ 
viding  the  spatial  and  temporal  resolution  necessary  to  char¬ 
acterize  a  variety  of  fundamental  issues  where  a  detailed  un¬ 
derstanding  is  still  lacking,  such  as  the  initial  chemical  steps 
of  decomposition,  detailed  chemical  reaction  pathways  (uni- 
and  multi-molecular),  and  the  molecular  structure  of  prod¬ 
ucts.  Such  fundamental  information  is  necessary  for  the  de¬ 
velopment  of  physics-based,  predictive  materials  models  that 
can  be  used  in  large-scale  macroscopic  simulations  of  HE 
materials. 

The  fundamental  input  in  MD  simulations  are  the  atomic 
interactions,  which  are  described  from  hrst  principles  by 
quantum  mechanics  (QM).  Unfortunately  QM  methods  are 
computationally  too  intensive  to  describe  most  processes  of 
interest  under  realistic  conditions  (especially  when  chemical 
reactions  need  to  be  described).  This  limitation  has  been 
overcome,  to  a  large  extent,  in  recent  years  with  the  devel- 
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opment  of  reactive  force  helds  (EEs)  that  can  describe 
chemical  reactions  in  a  computationally  efficient  way.^^^ 
ReaxEE"^  is  a  new-generation  bond  order-dependent  reactive 
EE  based  on  extensive  ab  initio  QM  calculations  that  enables 
for  the  hrst  time  the  accurate  simulation  of  HE  materials 
under  realistic  conditions  in  a  computationally  efficient  way. 

Typical  HE  materials,  such  as  TNT,  RDX,  HMX,  TATB, 
and  PETN,  are  complex  organic  molecular  solids;  their  deto¬ 
nation  involves  the  interplay  between  the  mechanical  shock 
loading,  the  induced  thermo-mechanical  response  (for  in¬ 
stance,  the  coupling  between  lattice  modes  and  internal  mo¬ 
lecular  modes,  large  internal  deformations  of  the  molecules, 
and  plastic  deformation)  and  the  induced  chemistry.  The  un¬ 
derstanding  of  detonation  was  advanced  signihcantly  with 
the  development  of  the  ZND  model  by  Zeldovich,  von  Neu¬ 
mann,  and  Doring  during  World  War  II.®  The  leading  shock 
front  compresses  and  heats  up  the  unreactive  material.  Eor 
strong  enough  shock  waves  the  high  temperatures  and  pres¬ 
sures  trigger  the  chemical  decomposition  of  the  HE  material 
and  is  followed  by  a  decrease  in  density  and  pressure  as  the 
chemical  reactions  progress  and  the  material  expands  in  the 
reaction  zone.®  Eor  steady  detonations  the  reaction  zone 
moves  with  the  leading  shock  at  the  detonation  velocity  and 
ends  in  the  so-called  Chapman-Jouget  (CJ)  point;  a  Taylor 
release  wave  follows,  where  the  products  expand  and  the 
pressure  decreases  to  a  value  that  depends  on  the  rear  bound¬ 
ary  condition  (equal  to  zero  for  unsupported  detonations).  A 
necessary  condition  for  the  steady  propagation  of  a  detona¬ 
tion  is  the  existence  of  a  “sonic  point”  that  isolates  the  re¬ 
action  zone  from  the  release  wave.  The  local  sound  speed 
decreases  from  right  behind  the  shock  front  into  the  reaction 
zone  and  Taylor  wave  due  to  expansion;  at  the  sonic  point 
the  local  sound  speed  is  equal  to  the  detonation  velocity  such 
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that  no  perturbation  in  the  Taylor  release  wave  can  affect  the 
reaction  zone.®  Real  HE  formulations  are  defective  (plastic 
bonded,  porous,  etc.)  and  these  heterogeneities  lead  to  hot¬ 
spots  and  a  complex  multi-dimensional  flow  even  in  macro- 
scopically  unidirectional  detonation. 

MD  with  reactive  potentials  that  exhibit  simple,  fast,  and 
exothermic  chemistry  has  been  extremely  useful  to  charac¬ 
terize  detonation  in  model  systems.  Nonequilibrium  MD 
shock  simulations  with  these  potentials  enabled  the  atomistic 
description  of  steady  detonation  waves, and  have  been 
used  to  study  the  role  of  defects  on  initiation,^’'®  the  exis¬ 
tence  of  a  failure  diameter,"  desensitization,*^  etc.  The  direct 
observation  of  steady  detonations  is  possible  because  these 
model  potentials  are  tuned  to  have  a  very  short  reaction  zone 
(around  10  run).  In  most  high  energy  materials  (including 
RDX,  HMX,  and  TATB)  the  reaction  zones  are  much  wider 
(up  to  millimeters)  and  reaction  times  much  longer  (a  frac¬ 
tion  of  a  microsecond),  making  direct,  nonequilibrium  atom¬ 
istic  simulations  of  steady  detonations  impossible  with  cur¬ 
rent  or  near-future  computational  capabilities.  Nevertheless 
reactive  MD  simulations  using  ReaxFF  provided  important 
information  regarding  the  initial  chemical  events  in  RDX 
under  shock  loading.®  NO2  was  found  to  be  the  only  product 
in  weak  shocks  in  agreement  with  experiments  and,  despite 
the  fact  that  the  barriers  for  pathways  leading  to  NO2  and 
HONO  are  essentially  the  same,  we  also  found  that  for 
strong  shocks  the  primary  reactions  leading  to  NO2,  OH, 
NO,  and  N2  occur  at  very  early  stages  with  time  scales  of  a 
few  ps. 

Equilibrium  MD  simulations  can  provide  valuable  infor¬ 
mation  regarding  the  thermal  decomposition  of  HE  materials 
enabling  the  characterization  of  processes  with  time  scales 
much  longer  than  those  attainable  in  nonequilibrium  shock- 
waves.  For  example,  Manaa  and  collaborators'®  used  QM- 
based  MD  to  study  HMX  decomposition  at  a  temperature 
and  density  (7’=  3500  K  and  p=1.9  g/cm®)  close  to  CJ  con¬ 
ditions.  Their  simulations  shed  light  into  the  initial  chemical 
reactions  and  time  scales  for  the  production  of  several  prod¬ 
ucts.  Unfortunately,  the  computational  cost  of  such  QM 
methods  hinders  their  applicability  to  study  decomposition  at 
lower  temperatures  that  require  much  longer  simulation 
times  (several  nanoseconds).  Another  promising  avenue  to 
explore  long  time-scale  phenomena  with  MD  is  the  use  of 
“hugoniostat”  techniques,  where  the  equations  of  motion  of 
the  atoms  are  modihed  to  reproduce  the  state  of  the  material 
caused  by  the  passage  of  a  shockwave. 

In  this  paper  we  use  ReaxFF  with  MD  to  study  thermal- 
induced  decomposition  of  RDX  at  various  densities  (from 
p=2.Il  g/cm®  in  compression  to  0.21  g/cm®  in  expansion) 
and  over  a  wide  temperature  range  (from  T=  1200  to  3000 
K).  Our  most  extreme  simulations  (p=2.11  g/cm®  and  T 
=  3000  K)  correspond  to  conditions  close  to  those  of  the  CJ 
point  (final  state  of  detonation)  of  RDX.  In  Sec.  II  we  briefly 
describe  ReaxFF  and  its  parametrization.  Section  III  de¬ 
scribes  our  MD  experiments  on  thermal  decomposition  and 
the  molecule  recognition  method  we  used  to  analyze  the 
simulations.  In  Secs.  IV-VI,  we  present  the  results  and  their 
implications.  Finally,  in  Sec.  VII  conclusions  are  drawn. 


II.  ReaxFF  FORCE  FIELD 

In  this  paper  we  used  the  ReaxFF  reactive  force  field  for 
nitramines  (ReaxFF„,,^o)  as  previously  used  to  study  the  ini¬ 
tial  chemical  events  in  RDX  during  shock  simulations.®  The 
potential  functions  in  ReaxFF,,, are  very  similar  to  those 
earlier  reported  for  the  hydrocarbon  reactive  force  field."' 
Here  we  will  only  describe  the  general  concepts  of  ReaxFF, 
followed  by  a  more  detailed  discussion  of  the  differences 
between  ReaxFF„,„.o  and  ReaxEE^/f . 

ReaxFF,,, partitions  the  system  energy  into  the  contri¬ 
butions  described  in  Eq.  (1): 

^  total  ^  bond~^  ^  over~^  ^  tinder'^  ^  lp~^  ^  triple~^  ^  val'^  ^  pen 

^  tors  ^  coa~^  ^  con]  ^vdW~^  ^  Co  u  ■  (^) 

A  fundamental  difference  between  the  ReaxFF  and  un¬ 
reactive  force  fields  is  that  ReaxFF  does  not  use  fixed  con¬ 
nectivity  assignments  for  the  chemical  bonds.  Instead  the 
bond  order,  BO',  is  calculated  directly  from  the  instanta¬ 
neous  interatomic  distances  r^j  [see  Eq.  (2)],  which  are  up¬ 
dated  continuously  during  the  dynamics.  This  allows  for  the 
creation  and  dissociation  of  bonds  during  a  simulation.  These 
instantaneous  bond  orders  BO'  are  corrected  for  overcoordi¬ 
nation  of  the  atoms  sharing  the  bond,  yielding  corrected 
bond  orders  BO  that  are  used  subsequently  on  all  covalent 
interactions  [e.g.,  the  functions  describing  bonds  'Ei,„„d, 
bond  angles  (E,,,,^,^  and  torsions  (E,o„)]: 

BO,^  =  BO,7+  BO,7+  BO,7’" 


/  \ «! 

exp 

p\ 

-l-exp 

P  7777 

7777 

l'‘o 

7777. 

All  covalent  interactions  are  expressed  in  terms  of  these 
bond  orders  so  that  these  terms  dissociate  properly  as  any 
bond  is  broken.  Since  bonds  can  break  and  form  during  dy¬ 
namics,  we  cannot  exclude  Coulomb  and  van  der  Waals  in¬ 
teractions  for  bonded  atoms  as  commonly  done  in  many  FF. 
Instead  we  must  include  these  interactions  between  every 
atom  pair,  independent  of  the  instantaneous  connectivity. 
This  is  very  natural  for  the  Coulomb  interactions  since  Re¬ 
axFF  treats  the  atomic  charges  as  having  the  finite  size  of  the 
atom,  leading  to  shielding  of  the  Coulomb  interaction  for 
shorter  distances.  We  also  include  a  similar  shielding  for  van 
der  Waals  interactions,  reducing  repulsion  at  short  inter¬ 
atomic  distances  (see  Eq.  12  in  Ref.  4). 

A.  New  terms  and  modifications  in  ReaxFF 

To  enable  the  description  of  nitramine  chemistry  the  fol¬ 
lowing  additions  were  made  to  ReaxFF^// ; 

Bond  energy  function.  Both  ReaxFF;-//  and 
ReaxFF„,,,.„  allow  separate  functional  forms  (BOg.,  BO,,, 
and  BO„.„)  to  describe  single,  double,  and  triple  bonds.  Each 
of  these  bond  orders  has  a  different  dependence  on  bond 
distance  with  the  parameters  determined  from  calculations 
on  molecules  such  as  H3C-NH3,  H2C=NH  and  HC=N. 
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However,  in  ReaxFFcn  the  total  bond  energy  is  described  as 
a  function  only  of  the  total  bond  order  (BOo.+  BO,,j 
+  BO^^).  In  contrast  ReaxFF„,-,„  uses  separate  energies  for 
single,  double,  and  triple  bonds,  as  shown  in  Eq.  (3).  The 
reason  for  this  change  is  that  while  carbon-carbon  systems 
have  bond  energies  that  increase  systematically  for  single  to 
double  to  triple  bonds  (dissociation  energies  of  98,  178,  235 
kcal/mol),  multiple  bonds  in  systems  containing  N  and  O  can 
have  dramatically  different  behaviors  as  bond  orders  in¬ 
crease.  This  arises  because  of  the  role  in  the  bonding  of  lone 
pairs  on  N  and  O,  which  can  participate  in  bonding.  Thus 
ReaxFF  has  been  extended  to  partition  the  total  bond  energy 
into  individual  contributions  from  sigma-bonds,  pi-bonds, 
and  double  pi-bonds  (triple  bonds),  as  indicated  in  Eq.  (3). 
As  in  ReaxEEc//,  the  bond  orders  BO'  as  obtained  from  Eq. 
(2)  are  corrected  for  overcoordination.  Using  these  corrected 
contributions  (BO),  Eq.  (3)  is  used  to  calculate  bond  energies 
in  ReaxEE„„,.o ; 


ij 


Lone  pair  electrons.  The  second  new  potential  function 
used  in  ReaxEE,,,,^^  is  an  energy  term  associated  with  lone 
electron  pairs  [Eip  in  Eq.  (1)].  Lone  pairs  of  electrons  gen¬ 
erally  play  little  role  in  hydrocarbon  chemistry,  but  lone  pairs 
on  such  heteroatoms  as  oxygen  and  nitrogen  are  important  in 
the  response  of  these  atoms  to  over-  and  undercoordination. 
Eurthermore,  the  presence  of  these  lone  electron  pairs  influ¬ 
ences  the  valence  angles  around  atoms  (this  aspect  was  al¬ 
ready  incorporated  in  ReaxFEc/^).  In  addition,  lone  pairs 
contribute  to  the  stability  of  conjugated  systems  by  delocal¬ 
izing  over  adjacent  pi  bonds.  Thus  Eq.  (4)  is  used  to  deter¬ 
mine  the  number  of  lone  pairs  around  an  atom.  A;  in  Eq.  (4) 
is  the  difference  between  the  total  number  of  outer  shell  elec¬ 
trons  (6  for  oxygen,  4  for  silicon,  1  for  hydrogen)  and  the 
sum  of  bond  orders  around  an  atomic  center. 


‘ip,i 


=  int 


[yj+exp 

^16 

2-f  Af-2int 


.  (4) 


Eor  oxygen  with  normal  coordination  (total  bond 
order=2,  A''  =  4),  Eq.  (4)  leads  to  two  lone  pairs.  As  the 
total  bond  order  associated  with  a  particular  O  atom  starts  to 
exceed  2,  Eq.  (4)  causes  a  lone  pair  to  gradually  break  up. 
This  is  accompanied  by  an  energy  penalty,  as  calculated  by 
Eq.  (5).  A  ip  i  in  Eq.  (5)  depicts  the  deviation  of  the  number 
of  lone  pairs  around  an  atom  from  the  number  of  lone  pairs 
at  normal  coordination  (2  for  oxygen,  1  for  nitrogen,  0  for 
carbon  and  hydrogen): 


,  _ _ Pip^ipj _ 

'P  \+exp{-15^,pp) 


(5) 


To  further  account  for  the  influence  of  lone  pairs  on 
molecular  stability  the  energy  expressions  used  in  ReaxEE^H 
to  calculate  the  contributions  of  overcoordination  energy 
have  been  expanded  to  account  for  deviations  in  the 
number  of  lone  pairs.  The  overcoordination  for  ReaxEE,,,,^^ 
differs  from  that  for  ReaxEE^H  in  that  the  overcoordination 
penalty  gets  diminished  for  systems  like  H2C=N=N  and 


R-N02  which  contain  an  over  coordinated  atom  as  a  result 
of  breaking  up  a  formal  lone  electron  pair  (A;^  ,  =  1)  next  to 
a  formally  undercoordinated  atom. 

Valence  angle  conjugation.  In  simple  valence  bond 
theory,  this  is  described  in  terms  of  a  formal  charge  transfer 
leading  to  an  N-f  atom  that  can  make  four  bonds  and  for 
H2C=N=N  an  N—  atom  that  can  make  only  two  bonds. 
ReaxEEc/f  contained  a  conjugation  term  keyed  to  torsion 
angles.  This  accounted  for  the  effect  of  conjugation  on  rota¬ 
tional  barriers.  However,  for  nitro  systems  the  main  effect  of 
N02-conjugation  is  on  the  stability  of  the  system  and  as  such 
cannot  describe  the  N02-conjugation.  To  properly  describe 
the  stability  of  the  nitro  group  we  now  include  a  conjugation 
term  keyed  to  valence  angles  (£^00)  to  the  ReaxEE  energy 
description.  Equation  (6)  shows  the  potential  function  used 
to  describe  the  valence  angle  conjugation  contribution  in 
angle  ijk: 

_  2  +  exp(-\2iAy) 

coa  ijk  i+exp(-\2iA2)  +  exp(\22A2) 

X  exp[  -  X  2o(  BOy  -  2 )  ]  exp[  -  X  2o(  BO^^t  -  2 )  ] , 

(6) 

where  A^  describes  the  overcoordination  on  the  central  atom 
in  valence  angle  ijk. 

Terminal  triple  bond  stabilization.  While  ReaxEE  rec¬ 
ognizes  triple  bonds,  a  bond  order  of  3.0  is  attained  at  dis¬ 
tances  shorter  than  equilibrium  since  it  reflects  an  asymptotic 
bond  order  value.  At  equilibrium,  a  system  with  a  formal 
triple  bond  usually  has  a  bond  order  of  about  2.7.  Although 
the  atoms  involved  in  the  triple  bond  do  not  completely  fill 
their  valency,  this  loss  of  bonding  is  countered  for  systems 
with  internal  triple  bonds  by  undercoodination  stabilization 
allowing  ReaxEE  to  properly  describe  C=C  bonds.  How¬ 
ever,  the  approach  used  in  ReaxEE^//  leads  to  too  little  sta¬ 
bility  in  cases  with  terminal  triple  bonds  as  allowed  by  N  and 
O  systems  (HC=N,  C=0,  and  N=N).  Eortunately,  it  is 
straightforward  to  detect  terminal  triple  bonds  allowing  us  to 
use  Eq.  (7)  to  account  for  the  additional  stabilization: 

E,np=v,rip  exp[-X8(BO,;-2.5)^] 

1 

^l-f25  exp[X5(A,  +  Ay)] 

X[exp[-X4(SBO,-BOj)] 

+  exp[-X4(SBO,.-BOj)]],  (7) 

SBO,  and  SBO^  are  the  sum  of  the  bond  orders  around  atoms 
i  and  j  and  A,  and  A^  describe  the  overcoordination  in  atoms 
i  and  j,  by  subtracting  the  valency  of  these  atoms  from  SBO, 
and  SBOy . 

All  the  modifications  introduced  in  ReaxEE^,,,.^  are 
fully  compatible  with  ReaxEE^//  and  with  the  ReaxEE  poten¬ 
tials  for  silicon/silicon  oxides*^  and  aluminum/aluminum 
oxides.*^ 

B.  Parametrization  of  ReaxFF 

The  parameters  of  the  nitramine  ReaxFF  are  based  on  a 
large  number  of  ab  initio  QM  calculations.  Over  40  reactions 
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FIG.  1.  Energetics  of  unimolecular  decomposition  mechanisms  in  RDX  obtained  using  the  ReaxFF  (full  lines  with  filled  symbols)  and  with  QM  (dashed  lines 
with  open  symbols)  (Ref.  20).  Circles  represent  the  sequential  HONO  elimination,  triangles  show  the  decomposition  process  following  homolytic  N-N  bond 
breaking  (NO2  elimination),  and  diamonds  represent  the  concerted  ring-opening  pathway.  Intermediates  and  products  are  described  in  Ref.  20. 


and  over  1600  equilibrated  molecules  have  been  used;  they 
are  designed  to  characterize  the  atomic  interactions  under 
various  environments  [likely  and  unlikely  (high  energy)] 
each  atom  can  encounter.  The  training  set  contains  bond 
breaking  and  compression  curves  for  all  possible  bonds, 
angle  and  torsion  bending  data  for  all  possible  cases,  as  well 
as  crystal  data.  The  data  used  to  parametrize  ReaxFF„,f^o 
includes  all  the  cases  previously  used  in  the  hydrocarbon 
training  set  described  in  Ref.  4  plus  the  additional  cases  de¬ 
tailed  in  the  supplementary  material  of  Ref.  5  which  also 
compares  the  accuracy  of  ReaxFF  in  fitting  the  QM  data. 

Thus  ReaxFF„,-„.„  accurately  describes  both  the  gas 
phase  chemistry  of  hydrocarbons  and  various  C,  H,  N,  O 
systems,  including  the  derived  decomposition  pathways  of 
RDX.  Figure  1  shows  ReaxFF  (full  lines)  and  QM  (dashed 
lines)  results  of  the  energetics  of  21  intermediate  and  transi¬ 
tion  states  involved  in  the  four  most  important  gas  phase 
decomposition  mechanisms.  We  show  the  main  three  mecha¬ 
nisms  studied  by  Chakraborty^'*  [sequential  HONO  elimina¬ 
tion,  homolytic  cleavage  of  an  NN  bond  (NO2  elimination) 
and  subsequent  decomposition,  and  concerted  decomposi¬ 
tion]  as  well  as  a  pathway  involving  NO2  insertion  suggested 
by  Zhang  et  al.  for  HMX.^^  For  this  last  case  we  performed 
QM  simulations  at  the  same  level  used  by  Chakraborty  et  al. 
for  the  RDX-equivalent  mechanism.  We  see  from  Fig.  1  that 
ReaxFF  accurately  describes  the  complex  gas  phase  decom¬ 
position  of  RDX. 


ReaxFF  also  describes  crystal  data  accurately;  we  obtain 
lattice  parameters  of  a  =  13.7781  A,  /7=12.0300A,  and  c 
=  10.9609  A  in  good  agreement  with  the  experimental  values 
(13.182,  11.574,  and  10.709  A).^^  The  calculated  bulk  modu¬ 
lus  Bj-=13.90GPa  is  also  in  good  agreement  with  the  ex¬ 
perimental  measurements  that  range  from  13.0  GPa^^  (iso¬ 
thermal)  to  11.1-11.4  GPa^"^  (isentropic,  from  ultrasonic 
resonant  frequencies).  Both  lattice  parameters  and  bulk 
modulus  were  obtained  by  fitting  Rose’s  equation  of  state^^ 
to  energy-volume  data  obtained  from  isothermal,  isocoric 
(NVT)  MD  at  7=  300  K.  The  lattice  parameters  for  the  NVT 
simulations  were  obtained  by  homogeneously  scaling  those 
resulting  from  energy  minimization  (relaxing  atomic  posi¬ 
tions  and  lattice  parameters). 

III.  MOLECULAR  DYNAMICS  SIMULATIONS 
OF  RDX  THERMAL  DECOMPOSITION 

We  study  the  decomposition  of  RDX  at  various  tempera¬ 
tures  between  T=  1200  K  and  r=3000K  and  at  three  densi¬ 
ties:  at  1.68  g/cm^  (close  to  normal  density — corresponding 
to  zero  pressure  volume:  V=Vq),  under  compression  (2.11 
g/cm^)  and  at  low  density  (0.21  g/cm^)  using  the  reactive 
potential  ReaxFF  with  MD.  The  initial  state  of  the  simula¬ 
tions  consist  of  RDX  perfect  crystals  using  simulation  cells 
containing  eight  molecules  (one  unit  cell,  168  atoms)  and 
3-D  periodic  conditions.  After  relaxing  the  atomic  positions 
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at  each  density  with  low  temperature  MD  (7’=10K)  we 
study  the  time  evolution  of  the  system  at  the  desired  tem¬ 
perature  with  isothermal  isochoric  (NVT  ensemble)  MD 
simulations  (using  a  Berendsen  thermostat;  the  time  scale 
associated  with  the  coupling  between  the  thermostat  and  the 
atomistic  system  was  set  to  200  fs).  In  order  to  correctly 
describe  the  fast  chemistry  under  extreme  conditions  we  use 
a  time  step  of  0.1  fs  in  integrating  the  equations  of  motion 
(nonreactive  simulations  can  typically  use  one  order  of  mag¬ 
nitude  larger  integration  step). 

Using  just  a  single  unit  cell  facilitated  the  exploration  of 
lower  temperatures  processes  where  simulations  for  several 
nanoseconds  were  necessary  to  reach  an  equilibrium  distri¬ 
bution  of  products  (one  nanosecond  corresponds  to  ten  mil¬ 
lion  MD  steps).  We  consider  that  use  of  one  unit  cell  (eight 
molecules)  is  sufficient  to  characterize  the  overall  time  evo¬ 
lution  of  the  chemistry  during  cook-off.  In  order  to  test  sys¬ 
tem  size  dependence  of  our  results  we  performed  some  simu¬ 
lations  with  larger  cells  (2X2X2  unit  cells,  64  RDX 
molecules,  and  1344  atoms).  Both  time  scales  and  equilib¬ 
rium  distribution  of  products  are  in  excellent  agreement  with 
the  smaller  simulations. 

Identifying  chemical  processes  under  the  extreme  condi¬ 
tions  characteristic  of  HE  decomposition  is  very  challenging 
due  to  the  high  energies  and  fast  time  scales  involved.  In 
order  to  discuss  these  chemical  processes  in  terms  of  mol¬ 
ecules  we  must  first  define  what  we  understand  to  be  a  mol¬ 
ecule  which  in  turn  is  based  on  the  definition  of  bond.  Bonds 
are  usually  dehned  in  configuration  space  (positions);  when 
two  atoms  are  closer  than  some  cutoff  distance  (depending 
on  the  elements  involved)  they  belong  to  the  same  molecule. 
However,  under  the  extreme  conditions  of  temperature  and 
pressure  studied  here  two  atoms  may  be  close  in  configura¬ 
tional  space  for  times  shorter  than  a  vibrational  period  (if 
their  c.m.  kinetic  energy  is  larger  than  the  binding  energy). 
Thus,  we  use  an  alternative  definition  that  two  atoms  are 
bonded  if  they  are  close  in  phase  space  (atomic  positions  and 
momenta);  in  practical  terms  we  require  the  two  atoms  to 
have  negative  relative  energy: 

Eu=yWirij)  +  K^i;!<0,  (8) 

where  is  the  effective  pairwise  interaction  between 

atoms  i  and  j  separated  by  a  distance  and  is  their 
c.m.  kinetic  energy.  The  interaction  between  two  atoms  in 
ReaxFF  depends  very  strongly  on  their  environment;  for  the 
purpose  of  defining  bonds  we  estimate  V\‘Jf(r)  using  an  ef¬ 
fective  two  body  interaction  with  a  modified  Morse  term: 


,,,,  \D^j{xHr)-2x(r))  if  r>R,j, 

\-Dij  if  r<Rij, 


(9) 


where  ;^'(r)  =  exp[7/2(l  —  r/R,^)].  In  order  to  define  unbi¬ 
ased  parameters  for  Eq.  (9)  the  energy  of  the  effective  inter¬ 
actions  {Djj)  are  proportional  to  the  ReaxFF  bond  energy 
parameter"^  and  the  distances  (Rjj)  are  proportional  to  the  cr 
bond  distances."^  The  three  free  parameters  [energy  scale,  dis¬ 
tance  scale,  and  curvature  (y)]  were  chosen  to  correctly  rec¬ 
ognize  RDX  molecules  in  their  crystalline  form  at  T 
=  300K.  The  prefactor  used  to  define  Djj  from  bond  ener¬ 


FIG.  2.  Time  evolution  of  potential  energy  (in  kcal/mol  per  unit  cell)  for 
temperatures  from  1200  to  T=3000  K  and  volumes  V=0.8Vo  (a),  V=Vo 
(b),  and  V=8Vo  (c). 


gies  is  5 ;  this  leads  roughly  to  the  energy  of  the  cr  bonds.^ 
The  distance  parameters  are  obtained  by  multiplying  the  cr 
bond  distances  by  1.35;  this  leads  to  equilibrium  bond  dis¬ 
tances.  In  order  to  avoid  having  different  molecules  in  the 
compressed  crystals  appear  to  be  clustered  together  we  chose 
a  large  value  for  the  curvature  7=25.  Similar  cluster  recog¬ 
nition  methods  have  been  used  in  fragmentation.^® 

IV.  ENERGETICS  AND  TIME  SCALES 
OF  DECOMPOSITION 

Figure  2  shows  the  time  evolution  of  the  potential  energy 
of  RDX  for  temperatures  T=  1200,  1500,  2000,  2500,  and 
3000  K  at  three  densities;  under  compression  (2.11  g/cm^) 
[Fig.  2(a)],  for  normal  density  (1.68  g/cm^)  [Fig.  2(b)],  and 
for  the  low  density  case  (0.21  g/cm^)  [Fig.  2(c)].  It  is  clear 
see  from  Fig.  2  that  higher  densities  lead  to  faster  chemistry. 
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TABLE  I.  Parameters  describing  the  exponential  behavior  of  potential  en¬ 
ergy  with  time  during  RDX  decomposition  obtained  from  our  MD  simula¬ 
tions  as  a  function  of  temperature  (T)  and  density  (p).  We  also  show  the 
equilibration-induction  time  for  each  case  Energies  are  given  in 

kcal/mol  and  times  in  ps. 


r(K) 

p  (g/cm^) 

tr-i 

Uo 

T 

1200 

2.11 

10 

-18  565.9 

2475.66 

1260.18 

1500 

2.11 

8 

-18  502.1 

2992.04 

248.28 

2000 

2.11 

3 

-18  043.6 

2541.12 

50.44 

2500 

2.11 

2.5 

-17  622.8 

3160.93 

14.47 

3000 

2.11 

2 

-17  355.0 

3217.71 

6.94 

1200 

1.68 

10 

-18  667.6 

2610.23 

1857.16 

1500 

1.68 

9 

-18  423.3 

2759.97 

276.96 

2000 

1.68 

5 

-18213.3 

3150.71 

59.96 

2500 

1.68 

2 

-17931.4 

3218.99 

17.73 

3000 

1.68 

1 

-17  720.0 

3201.26 

10.73 

1500 

0.21 

25 

-17  705.5 

2102.23 

2688.19 

2000 

0.21 

22 

-17  004.9 

2032.99 

498.93 

2500 

0.21 

17.5 

-16  929.8 

2494.95 

169.34 

3000 

0.21 

15.0 

-16  743.3 

2347.87 

112.44 

After  an  initial  equilibration  and  induction  time  (denoted 
?£_/)  the  potential  energy  decreases  with  time  as  the  chemi¬ 
cal  reactions  progress.  Its  temporal  evolution  can  then  be 
described  rather  well  with  a  simple  exponential  function; 

Uif,T,p)  =  Uo{T,p}  +  AU(T,p) 

Xexp[-(t-tE_j)/T{T,p)],  (10) 

where  UqIS  the  asymptotic  energy  of  the  products,  A  t/  is  the 
exothermicity  of  the  reaction,  and  t  is  an  overall  character¬ 
istic  time  of  reaction;  all  three  parameters  are  temperature 
and  density  dependent.  In  Table  I  we  show  the  parameters 
obtained  from  fitting  Eq.  (10)  to  our  MD  data. 

Figure  3  shows  the  characteristic  times  (in  logarithmic 
scale)  as  a  function  of  inverse  temperature  (Arrhenius  plot). 
The  dashed  line  corresponds  to  MD  results  at  2.11  g/cm^,  the 
solid  line  with  open  squares  corresponds  to  1.68  g/cm^,  and 
the  dotted  one  represents  the  theoretical  results  0.21  g/cm^. 
The  solid  line  shows  the  behavior  of  HMX  (an  HE  material 


FIG.  3.  Characteristic  time  versus  inverse  temperature  from  ReaxFF  MD 
simulations  for  different  volumes:  V=0.8Vo  (red  line),  V=Vo  (black  line), 
and  V=8Vo  (green  line).  We  also  include  the  Arrhenius  behavior  obtained 
from  a  wide  range  of  experimental  ignition  times  in  HMX  (blue  line) 
(Ref.  27). 


very  similar  to  RDX)  obtained  from  a  wide  variety  of  experi¬ 
mental  ignition  times,  including  fast  processes  (10^®  s),  by 
Henson  and  collaborators.^^  The  compilation  of  experimental 
data  includes  a  wide  variety  of  sources  such  as  thermal  ex¬ 
plosions,  fast  pyrolysis,  laser  heating,  and  detonations.^^  Fig¬ 
ure  3  shows  that  our  first  principles-based  calculations  are  in 
reasonable  quantitative  agreement  with  experiments.  The 
Arrhenius  behavior  shown  by  both  MD  and  experimental 
results  can  be  attributed  to  a  single  rate-limiting  step  in  the 
sequence  of  reactions  that  lead  to  decomposition.  From  the 
MD  simulations  we  obtain  an  activation  energy  {E^)  for  this 
limiting  step  that  increases  as  the  density  decreases:  for 
p=2.11  g/cm^  we  obtain  £^  —  23  kcal/mol  and  for  p=0.21 
g/cm^  £„  =  26.6  kcal/mol. 

V.  TIME  EVOLUTION  OF  INTERMEDIATES 
AND  PRODUCTS 

In  Fig.  4  we  show  the  time  evolution  of  the  population  of 
key  products  N2 ,  H2O,  CO,  and  CO2  and  intermediate  NO2 ; 
the  time  is  scaled  with  the  corresponding  characteristic  reac¬ 
tion  time  obtained  from  the  energy  evolution.  We  find  that 
N-N  bond  breaking  to  form  NO2  is  the  dominant  early 
chemical  reaction  for  the  conditions  studied.  This  occurs  de¬ 
spite  the  similar  energy  barrier  for  HONO  elimination.  This 
difference  is  expected  since  the  NO2  formation  mechanism 
leads  to  greatly  increased  entropy  (and  hence  less  free  en¬ 
ergy),  while  HONO  elimination  has  a  tight  transition  state 
(low  entropy  and  high  free  energy).  This  was  also  observed 
in  the  pressure  and  temperature  dependence  in  the  reaction 
rates  estimated  in  Ref.  28.  It  was  also  apparent  in  the  shock 
induced  decomposition^  whereas  NO2  products  dominate  the 
early  species  for  weak  shocks  (up  to  particle  velocity  of  3 
km/s)  while  higher  velocities  led  to  equal  amounts  of  HO 
and  NO  at  early  times.  For  every  condition  of  temperature 
and  density  we  find  that  the  time  scales  of  H2O  and  N2 
formation  are  similar  to  one  another  and  to  the  overall  char¬ 
acteristic  time  (r);  formation  of  H2O  and  N2  occurs  during 
the  early  stages  of  the  process.  Manaa  et  al.  found  that  H2O 
forms  faster  than  N2  in  HMX  decomposition  using  ab  initio 
MD.  We  also  find  that  the  time  scales  of  CO  and  CO2  for¬ 
mation  depend  on  density;  for  normal  and  high  density  the 
time  scales  of  CO  and  CO2  formation  are  larger  than  those  of 
H2O  and  N2  whereas  at  low  density  CO  forms  promptly  and 
with  time  scales  similar  to  those  of  H2O  and  N2  (see  Fig.  4). 

VI.  EQUILIBRIUM  STATE  OF  THE  PRODUCTS 
A.  Asymptotic  population  of  key  moiecuies 

Figure  5  shows  the  asymptotic  (time  averaged)  popula¬ 
tions  of  the  key  products  (N2,  H2O,  CO,  and  CO2)  for  the 
various  temperatures  and  densities  studied.  We  find  that  the 
population  of  N2  decreases  with  compression.  At  low  density 
essentially  all  N  atoms  form  N2  molecules;  for  p= 2. 11  g/cm^ 
the  population  of  N2  is  smaller  (involving  roughly  half  the  N 
atoms)  and  increases  with  temperature.  The  population  of 
H2O  also  decreases  under  compression  but  the  values  for 
p=0.21  and  1.68  g/cm^  are  comparable.  The  populations  of 
CO  and  CO2  show  a  much  more  marked  density  dependence. 
We  see  from  Fig.  5  that  at  low  densities  (p=0.21  g/cm^) 
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FIG.  4.  (Color)  Time  evolution  of  products  Ni ,  H2O,  CO,  and  CO2  and  intermediate  NO2  for  T=3000  K  and  T=  1500  K  and  three  densities  p=2. 1 1  g/cm^, 
p=1.68  g/cm^,  and  p=0.21  g/crc?.  Here  (b)  corresponds  approximately  to  the  CJ  point. 


virtually  all  C  atoms  form  CO  molecules  while  we  find  no 
CO2  molecules  (the  total  number  of  C  atoms  in  the  system  is 
24).  Figure  6(a)  shows  a  snapshot  of  the  equilibrium  state  for 
r=3000K  and  p=0.21  g/cm^.  All  atoms  form  small  mol¬ 
ecules;  in  addition  to  the  ones  shown  in  Fig.  5  we  find  O2 
and  H2  with  asymptotic  populations  of  4  and  7,  respectively. 
As  the  density  increases  the  equilibrium  population  of  CO2 
increases  at  the  expense  of  CO  molecules;  this  happens  be¬ 
cause  at  higher  densities  C  atoms  tend  to  clusterize  into  con¬ 
densed  phase  aggregates,  leading  to  a  C-poor  gas  phase.  This 
results  in  an  increase  in  the  population  of  CO2  molecules. 
For  p=2.11  g/cm^  we  find  almost  no  CO  molecules;  all  the  C 
atoms  that  do  not  form  CO2  molecules  end  up  in  the  con¬ 
densed  phase  aggregate.  Figure  6(b)  shows  a  snapshot  of  an 


equilibrium  configuration  for  y=0.8yo,  containing  a  small 
condensed  phase  cluster  with  ten  carbon  atoms. 

The  phenomena  of  carbon  clustering  is  very  important  to 
understanding  the  properties  of  HE  materials.  In  Fig.  7  we 
plot  the  maximum  number  of  C  atoms  in  a  single  molecule 
as  a  function  of  temperature  and  for  the  three  densities  stud¬ 
ied  here.  For  p=0.21  glcvc?  there  are  no  molecules  with 
more  that  one  C  atom  at  any  temperature.  For  p=1.68  g/cm^ 
and  high  temperatures  the  largest  carbon  aggregate  slowly 
grows  with  decreasing  temperature,  but  between  T=  1500 
and  1200  K  we  find  an  abrupt  jump  in  the  size  of  the  con¬ 
densed  phase.  For  p=2.11  g/cm^  we  find  larger  C  clusters  at 
all  temperatures  and  the  abrupt  increase  in  the  size  of  the  C 
aggregate  occurs  at  higher  temperatures  (between  2000  and 
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FIG.  6.  (Color)  (a)  Snapshot  of  an  equilibrium  configuration  for  T=3000  K 
and  p=0.21  g/cm^’.  (b)  Snapshot  for  T=3000  K  and  p=2. 11  g/cm^;  this 
coiTesponds  approximately  to  the  CJ  conditions,  showing  molecular  charac¬ 
ter  at  this  point. 


FIG.  5.  Final  population  of  key  molecules  for  various  temperatures  from 
T=3000  K  to  T=1200  K  and  three  densities:  p=2.11  g/cm^  (a),  p=1.68 
g/cm^  (b),  and  p=0.21  g/cm^  (c). 


2500  K).  Larger  systems  and  longer  simulations  are  neces¬ 
sary  to  fully  characterize  the  relatively  slow  process  of  car¬ 
bon  clustering,  nevertheless  our  simulations  seem  to  indicate 
the  existence  of  two  distinct  regimes:  one  dominated  by  large 
carbon  aggregates  (favored  by  high  density  and  low  tempera¬ 
ture)  and  a  second  in  which  most  C  atoms  form  small  mol¬ 
ecules  (for  low  densities  and  high  temperatures). 

B.  Mean  lifetimes  of  products  in  equilibrium 

Figures  4  and  5  show  total  populations  of  various  mol¬ 
ecules  but  do  not  provide  any  information  regarding  their 
atomic  stability:  what  is  the  average  lifetime  of  an  individual 
molecule?  Under  these  extreme  conditions  we  expect  popu¬ 
lations  to  be  in  dynamical  equilibrium  by  undergoing  fre¬ 
quent  atomic  rearrangements. 

A  direct  calculation  of  the  mean  lifetime  of  a  given  mol¬ 
ecule  type  would  involve  averaging  the  time  elapsed  be¬ 
tween  the  formation  of  a  given  molecule  and  its  disappear¬ 


ance.  This  direct  method  has  two  drawbacks:  First,  it  leads  to 
very  inaccurate  results  when  the  mean  lifetime  is  longer  than 
or  comparable  to  the  simulation  time.  The  second  problem  is 
that  it  is  difficult  to  recognize  unambiguously  molecules  in  a 


FIG.  7.  Maximum  number  of  C  atoms  in  a  single  molecule  as  a  function  of 
temperature  for  p=2.11  g/cm^  (squares),  1.68  g/cm^  (circles),  and  0.21 
g/cm^  (diamonds). 
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single  snapshot  in  time  of  atomic  positions  and  velocities.  A 
molecule  may  seem  to  disappear  (large  vibrational  ampli¬ 
tude)  and  reappear  shortly  after  (within  one  or  two  vibra¬ 
tional  periods).  Is  one  to  count  such  an  event  as  a  real  chemi¬ 
cal  reaction  or  should  it  be  considered  an  artifact  and 
discarded?  The  molecule  recognition  algorithms  used  in  MD 
have  limitations  in  properly  accounting  for  such  effects, 
leading  to  serious  problems  in  calculating  mean  lifetimes. 

We  use  a  more  robust  method  to  estimate  the  mean  life¬ 
time  of  molecules  from  MD  simulations  based  on  the  fact 
that  in  most  cases  several  statistically  equivalent  molecules 
of  the  same  type  are  present  at  any  given  time.  We  find  it 
convenient  to  evaluate  the  mean  lifetimes  ( T^)  of  molecule 
a  (once  its  population  reached  its  asymptotic  values)  using 


T  = 

^  rv 


Un'^-nS 


(11) 


where  At  is  the  period  of  time  used  in  the  calculation,  N „  is 
the  average  number  of  molecules  of  type  a,  and  N'^  is  the 
number  of  different  (atomistically  distinguishable)  molecules 
of  type  a  that  were  present  during  the  time  of  observation. 
The  numerator  in  Eq.  (11)  represents  to  total  effective  obser¬ 
vation  time;  the  denominator  represents  the  number  of  times 
a  molecule  of  type  a  is  created  or  destroyed.  If  N'^=N no 
reactions  occurred  during  the  observation  time  and  our  best 
estimate  of  the  mean  lifetime  is  infinity.  As  another  example 
of  the  meaning  of  Eq.  (11),  suppose  we  start  the  calculation 
with  ten  molecules  of  a  given  type  and  one  of  them  disap¬ 
pears  while  a  new  one  is  formed  during  the  observation  time 
(Af);  the  average  number  of  molecules  (A)  is  ten  and  the 
number  of  distinguishable  molecules  is  11.  The  mean  life¬ 
time  according  to  Eq.  (11)  is  r=10Af/2.  This  is  the  result 
we  would  obtain  if  we  observed  a  single  molecule  for  a  time 
10 At  and  our  molecule  disappears  and  a  new  one  is  formed 
during  the  observation  time.  Note  that  Eq.  (11)  is  only  valid 
once  the  molecular  population  has  reached  its  asymptotic 
value;  if  more  molecules  of  type  a  are  being  formed  than 
destroyed,  Eq.  (11)  will  underestimate  their  mean  lifetime. 

Eigure  8  shows  the  mean  free  lifetimes  for  N2,  H2O, 
CO,  and  CO2  as  a  function  of  inverse  temperature  for  the 
three  densities  studied.  Note  that  the  populations  of  CO  at 
p=2.11  g/cm^  and  that  of  CO2  at  p=0.21  g/cm^  are  very 
small  and  the  calculation  of  mean  lifetime  becomes  very  in¬ 
accurate.  Even  for  the  other  cases  the  statistical  uncertainties 
(due  to  small  system  size  and  short  time-scales  characteristic 
of  MD  simulation)  are  significant.  Our  MD  simulations  show 
that  for  all  conditions  studied  we  have  well  defined  mol¬ 
ecules  with  mean  lifetimes  between  —10  ps  and  10  ns.  The 
simulation  at  p=2.11  g/cm^  and  7=  3000  correspond  to  con¬ 
ditions  close  to  those  of  the  C-J  point  in  RDX.  We  find  that 
under  those  conditions  the  mean  lifetime  of  H2O  molecules 
is  —5  ps,  that  of  CO2  is  around  30  ps,  and  N2  lives  much 
longer,  —450  ps.  This  indicates  that  at  the  C-J  point  the 
system  is  in  a  molecular  state  (rather  than  a  plasma  like  state 
with  no  covalent  bonding). 


FIG.  8.  Equilibrium  mean  life  times  of  products  [Eq.  (11)]  as  a  function  of 
temperature  for  the  three  densities  studied. 


VII.  CONCLUSIONS 

In  summary,  we  studied  the  thermal  decomposition  of 
RDX  for  a  variety  of  temperatures  and  densities  using  the 
first  principles-based  reactive  force  field  ReaxEE  with  MD. 
Erom  our  simulations  we  extract  a  characteristic  reaction 
time  that  shows  Arrhenius  temperature  dependence  and  de¬ 
creases  with  density.  This  suggests  that  a  single  step  kinetics 
is  a  good  approximation  in  this  case.  We  also  studied  the 
time  evolution  and  asymptotic  population  of  various  prod¬ 
ucts  and  intermediates.  The  population  of  CO  and  CO2  show 
a  marked  density  dependence.  Eor  low  density  essentially  all 
C  atoms  form  CO  molecules,  with  an  asymptotic  population 
of  CO2  of  almost  zero.  With  compression  C  atoms  tend  to 
clusterize  into  a  condensed  phase  aggregate  leading  to  a 
C-poor  gas  phase  and  an  increase  in  the  population  of  CO2  at 
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the  expense  of  CO.  For  p=2.11  g/cm^  we  find  almost  no  CO 
molecules.  Our  simulations  shed  some  light  into  the  initial 
stages  of  the  process  of  carbon  clustering. 

The  molecular  state  of  the  material  at  CJ  conditions  is  a 
critical  question  in  the  energetic  materials  community  and 
very  hard  to  answer  experimentally.  Our  results  show  that 
under  conditions  near  to  the  CJ  point  (p=2.11  g/cm^  and  T 
=  3000  K)  molecules  with  relatively  long  mean  lifetimes  (in 
terms  bond  vibrations)  are  clearly  identifiable:  the  mean  life¬ 
time  of  N2  molecules  is  around  300  ps,  while  H2O  and  CO2 
have  lifetimes  of  5  and  30  ps,  respectively.  This  shows  that 
under  CJ  conditions  RDX  corresponds  to  a  molecular  state  as 
opposed  to  being  plasmalike  with  no  molecular  entities. 

The  accurate  description  of  the  shock  velocity  as  a  func¬ 
tion  of  piston  velocity  in  RDX  shocks^  and  the  reasonable 
time  scales  for  the  thermal  decomposition  of  RDX  (as  com¬ 
pared  with  experiments  in  a  similar  energetic  material, 
HMX)  indicate  that  ReaxFF  yields  an  accurate  description  of 
the  atomic  interactions  of  RDX  even  under  the  extreme  con¬ 
ditions  attained  during  detonation. 

Atomistic  modeling  with  these  new  generation  reactive 
potentials  is  becoming  an  increasingly  important  tool  to  in¬ 
vestigate  condensed-phase  chemistry  under  dynamical  or 
static  loading.  Such  simulations  provide  very  detailed  infor¬ 
mation  regarding  the  decomposition  and  subsequent  reac¬ 
tions  in  energetic  materials  that  allow  one  to  make  some 
sense  out  of  the  very  complex  mechanical  and  chemical  pro¬ 
cesses.  Furthermore,  since  ReaxFF  is  based  on  ab  initio  cal¬ 
culations,  all  the  predictions  in  this  paper  can  be  considered 
as  first  principles  (no  critical  empirical  data  was  used  in  our 
calculations).  Finally,  the  data  presented  in  this  paper  should 
provide  valuable  input  to  mesoscopic  and  macroscopic  mod¬ 
els  of  energetic  materials  and  these  type  of  simulations  are  an 
important  step  to  develop  physics-based,  predictive  materials 
models. 


ACKNOWLEDGMENTS 

Work  at  Los  Alamos  National  Laboratory  was  supported 
by  the  ASC  Materials  and  Physics  Modeling  Program, 
LANL.  Work  at  Caltech  was  supported  by  ONR  (program 
manager  Judah  Goldwasser). 


‘D.  D.  Dlott,  Annu.  Rev.  Phys.  Chem.  50,  251  (1999). 

^S.  D.  McGrane,  D.  S.  Moore,  and  D.  J.  Funk,  J.  Appl.  Phys.  93,  5063 
(2003). 

^D.  W.  Brenner,  0.  A.  Shenderova,  J.  A.  Harrison,  S.  J.  Stuart,  B.  Ni,  and 
S.  B.  Sinnott,  J.  Phys.:  Condens.  Matter  14,  783  (2002);  S.  J.  Stuart,  A.  B. 
Tutein,  and  J.  A.  Hanison,  J.  Chem.  Phys.  112,  6472  (2000). 

'^A.  C.  T.  van  Duin,  S.  Dasgupta,  F.  Lorant,  and  W.  A.  Goddard  III,  J.  Phys. 
Chem.  A  105,  9396  (2001). 

^A.  Strachan,  A.  C.  T.  van  Duin,  D.  Chakraborty,  S.  Dasgupta,  and  W.  A. 
Goddard  III,  Phys.  Rev.  Lett.  91,  098301  (2003). 

*W.  Picket  and  W.  C.  Davis,  Detonation  (Univ.  of  California,  Berkeley, 
1979). 

’M.  L.  Elert,  D.  M.  Deaven,  D.  W.  Brenner,  and  C.  T.  White,  Phys.  Rev.  B 
39,  1453  (1989). 

*D.  W.  Brenner,  D.  H.  Robertson,  M.  L.  Elert,  and  C.  T.  White,  Phys.  Rev. 
Lett.  70,  2174  (1993). 

'^J.  W.  Mintmire,  J.  J.  C.  Ban'ett,  D.  H.  Robertson,  and  C.  T.  White,  Chem. 
Phys.  Rep.  17,  37  (1998). 

'"B.  L.  Holian,  T.  C.  Germann,  J.-B.  Maillet,  and  C.  T.  White,  Phys.  Rev. 
Lett.  89,  285501  (2002). 

"C.  T.  White,  D.  H.  Robertson,  D.  R.  Swanson,  and  M.  L.  Elert,  AIP  Conf. 
Proc.  505,  377  (2000). 

'^B.  M.  Rice,  W.  Mattson,  and  S.  F.  Trevino,  Phys.  Rev.  E  57,  5106  (1998). 
'^M.  R.  Manaa,  L.  E.  Fried,  C.  F.  Melius,  M.  Elstner,  and  Th.  Frauenheim, 
I.  Phys.  Chem.  A  106,  9024  (2002). 

*'*1.-B.  Maillet,  M.  Mareschal,  L.  Soulard,  R.  Ravelo,  P.  S.  Lomdahl,  T.  C. 

Germann,  and  B.  L.  Holian,  Phys.  Rev.  E  63,  016121  (2001). 

*^E.  J.  Reed,  L.  E.  Fried,  and  J.  D.  Joannopoulos,  Phys.  Rev.  Lett.  90, 
235503  (2003). 

“•I.  K.  Brennan  and  B.  M.  Rice,  Mol.  Phys.  101,  3309  (2003). 

'^R.  Ravelo,  B.  L.  Holian,  T.  C.  Germann,  and  P.  S.  Lomdahl,  Phys.  Rev.  B 
70,  014103  (2004). 

'*A.  C.  T.  van  Duin,  A.  Strachan,  S.  Stewman,  Q.  S.  Zhang,  X.  Xu,  and  W. 

A.  Goddard  III,  J.  Phys.  Chem.  107,  3803  (2003). 

'®Q.  Zhang,  T.  Cagin,  A.  C.  T.  van  Duin,  W.  A.  Goddard  III,  Y.  Qi,  and  L.  G. 

Hector,  Phys.  Rev.  B  69,  045423  (2004). 

^°D.  Chakraborty,  R.  P.  Muller,  S.  Dasgupta,  and  W.  A.  Goddard  III,  I.  Phys. 
Chem.  A  104,  2261  (2000). 

^'S.  Zhang,  H.  N.  Nguyen,  and  T.  N.  Truong,  J.  Phys.  Chem.  A  107,  2981 
(2003). 

^“C.  S.  Choi  and  E.  Prince,  Acta  Crystallogr.,  Sect.  B:  Struct.  Crystallogr. 
Ciyst.  Chem.  B28,  2857  (1972). 

^^B.  Olinger,  B.  Roof,  and  H.  Cady,  Symposium  International  Sur  Le  Com¬ 
portment  Des  Milieux  Denses  Sous  Hautes  Pressions  Dynamiques,  Com¬ 
missariat  a  TEnergie  Atomique  Centre  d’Etudes  de  Vajours,  Paris,  France, 
1978,  p.  3. 

-'‘S.  Haussul,  Z.  Kristallogr.  216,  339  (2001). 

^^1.  H.  Rose,  J.  Ferrante,  and  I.  R.  Smith,  Phys.  Rev.  Lett.  47,  675  (1981). 
^®A.  Strachan  and  C.  O.  Dorso,  Phys.  Rev.  C  56,  995  (1997);  S.  Pratt,  C. 
Montoya,  and  F.  Romming,  Phys.  Lett.  B  349,  261  (1995);  J.  Pan  and  S. 
D.  Gupta,  Phys.  Rev.  C  51.  1384  (1995). 

^’B.  F.  Henson,  L.  Smilowitz,  B.  W.  Asay,  P.  M.  Dickson,  and  P.  M.  Howe, 
Proceedings  of  the  12th  Symposium  (International)  on  Detonation  (2002), 
http://www.sainc.com/detsymp/technicalProgram.htm 
“*D.  Chakraborty,  R.  P.  Muller,  S.  Dasgupta,  and  W.  A.  Goddard  III,  J. 
Comput. -Aided  Mater.  Des.  8,  203  (2001). 


Downloaded  12  Mar  2005  to  131.215.16.37.  Redistribution  subject  to  AiP  iicense  or  copyright,  see  http://jcp.aip.org/jcp/copyright.jsp 


